// morelos_turbo.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"


// -----------------------------------------------------------------------
// File:    turbo.cpp
// Date:    May 5, 2002
// Authors: M. Walma and R.H. Morelos-Zaragoza
//
// Description: Simulation of a binary rate-1/3 turbo code with rate-1/2
// component recursive convolutional codes.
//
// Based on a program written by Mathys Walma, who placed it in the public 
// domain in 1998. The user interface was modified by R. Morelos-Zaragoza.
// ------------------------------------------------------------------------
// This program is complementary material for the book:
//
// R.H. Morelos-Zaragoza, The Art of Error Correcting Coding, Wiley, 2002.
//
// ISBN 0471 49581 6
//
// This and other programs are available at http://the-art-of-ecc.com
//
// You may use this program for academic and personal purposes only. 
// If this program is used to perform simulations whose results are 
// published in a journal or book, please refer to the book above.
//
// The use of this program in a commercial product requires explicit
// written permission from the author. The author is not responsible or 
// liable for damage or loss that may be caused by the use of this program. 
//
// Copyright (c) 2002. Robert H. Morelos-Zaragoza. All rights reserved.
// ------------------------------------------------------------------------
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <limits.h>
#include "rng.h"


#define MAXRAND       LONG_MAX

int          states;
long         seed;
unsigned int N;
unsigned int ITERATIONS;
long         maxsim;
double       No;
double       sigma;
int          temp;

double *L[2];


double gaussian(double variance);
void encode(bool *mesg,bool *parity1,unsigned size, bool force);
void interleave(bool *mesg, unsigned size);
bool add(bool a,bool b);
void interleave(bool *mesg, unsigned size);
void deinterleave(bool *mesg, unsigned size);
void deinterleavedouble(double *mesg, unsigned size);
void interleavedouble(double *mesg, unsigned size);
void createencodetable(int states);
void inttobool(unsigned state, bool *array, unsigned size);
void booltoint(bool *array, unsigned size, unsigned *state);
unsigned decode (double *channelmesg, double *parity1, double *parity2, unsigned size, bool *mesg);

unsigned *interleavearray;
unsigned *deinterleavearray;

unsigned numstates;
unsigned memory;

// [2] = input, [16] = current state, tostate[2][16] = next state
unsigned *tostate[2];

// [2] = last input, [16] = current state, fromstate[2][16] = previous state
unsigned *fromstate[2];

// [2] = input, [16] = current state, output[2][16] = output of encoder
bool *output[2];

RNG gRandom;

bool boolrandom(){
    double num = gRandom.rand_closed01() * 2;
    int coin = (int)num;
    return coin ? true : false;
}

int _tmain(int argc, char* argv[])
{
   bool *mesg;
   bool *parity1;
   bool *parity2;
   double *channelmesg;
   double *channelparity1;
   double *channelparity2;
   register int i;
   char name1[40], name2[40];
   FILE *fp;
   gRandom.init(40539119,40639968,33552418,96088923);


   // Command line processing
   //if (argc != 10)
   //  {
   //    printf("Usage: %s Num_states Max_iter INT_file", argv[0]);
   //    printf(" BER_file EbNo_init EbNo_final EbNo_inc  max_sims  seed\n");
   //    exit(0);
   //  }

   double INIT, FINAL, INC;

   /*sscanf(argv[1],"%d",  &states);
   sscanf(argv[2],"%d",  &ITERATIONS);
   sscanf(argv[3],"%s",  name1);
   sscanf(argv[4],"%s",  name2);
   sscanf(argv[5],"%lf", &INIT);
   sscanf(argv[6],"%lf", &FINAL);
   sscanf(argv[7],"%lf", &INC);
   sscanf(argv[8],"%ld", &maxsim);
   sscanf(argv[9],"%ld", &seed);*/


   states=16;
   ITERATIONS=8;
   strcpy(name1,"int1.txt");;
   strcpy(name2,"ber.txt");
   INIT = 0.0;
   FINAL = 0.7;
   INC = 0.1;
   maxsim = 100000000;
   seed = 0xBABACA;
   _set_SSE2_enable(1);


   // READ INTERLEAVER
   fp = fopen(name1,"r");
   N = 0;
   while (!(fscanf(fp,"%d", &temp)==EOF))
     N++;
   fclose(fp);
   interleavearray = new unsigned[N];
   deinterleavearray = new unsigned[N];
   fp = fopen(name1,"r");
   i = 0;
   while (!(fscanf(fp,"%d", &interleavearray[i])==EOF))
     i++;
   fclose(fp);
   for (i=0; i<N; i++) deinterleavearray[interleavearray[i]]=i;

   double rate = (double)N / (3.0*(double)N + 4.0);

   printf("\nParallel concatenated (turbo) code of rate = %lf\n", rate);
   printf("%2d-state constituent recursive convolutional codes \n\n", states);

   srand(seed);

   mesg = new bool[N];
   parity1 = new bool[N];
   parity2 = new bool[N];
   channelmesg = new double[N];
   channelparity1 = new double[N];
   channelparity2 = new double[N];

   fp = fopen(name2,"w");

   FILE* sigs = fopen("sigs.txt","rb");
   FILE* signals = fopen("sigs_out.txt","wb");

   createencodetable(states);



   for (double dB=INIT; dB<FINAL+1.0e-1; dB += INC)
   {
	unsigned totalN=0;
	unsigned totalerrors=0;
	unsigned numiter;
	unsigned totaliter=0;
	unsigned numloops=0;
    unsigned totalglitches=0;
    unsigned totalbetters=0;

      	// dB = 10*log(Eb/No) where Eb is 1
      	// No = 1.0/pow(10.0,dB/10.0);
      	sigma = 1.0/sqrt(2.0*0.3333333*pow(10.0,dB/10.0));
	
      	bool keepgoing=true;

	do {
		for (int x=0;x<N;x++)
			mesg[x] = boolrandom();


         
       	// create a random interleaver for each decode trial
		// createinterleave(N);

		encode(mesg,parity1,N,true);
		interleave(mesg,N);
		encode(mesg,parity2,N,false);
		deinterleave(mesg,N);

        //fread(mesg,sizeof(bool),N,sigs);
        //fwrite(mesg,sizeof(bool),N,signals);


		for (int x=0;x<N;x++)
		   {
		   mesg[x]    ?    channelmesg[x]=1.0 :    channelmesg[x]=-1.0;
		   parity1[x] ? channelparity1[x]=1.0 : channelparity1[x]=-1.0;
		   parity2[x] ? channelparity2[x]=1.0 : channelparity2[x]=-1.0;
		   }

        
		// add gaussian noise, mean=0, variance=No/2
		for (int x=0;x<N;x++)
		   {
		    	//channelmesg[x] += gaussian(No/2);
			    //channelparity1[x] += gaussian(No/2);
		      	//channelparity2[x] += gaussian(No/2);
                /*channelmesg[x] += gRandom.normal(1,sigma);
			    channelparity1[x] += gRandom.normal(1,sigma);
		      	channelparity2[x] += gRandom.normal(1,sigma);*/
               channelmesg[x] += gaussian(sigma);
               channelparity1[x] += gaussian(sigma);
               channelparity2[x] += gaussian(sigma);
		   }

     /*   double power=0;
        for (int x=0;x<N;x++)
		{
               power += channelmesg[x] * channelmesg[x];
        }
        power /= N;
        double rms = sqrt(power);*/

        /*fread(channelmesg,sizeof(double),N,sigs);
        fread(channelparity1,sizeof(double),N,sigs);
        fread(channelparity2,sizeof(double),N,sigs);*/

        //fwrite(channelmesg,sizeof(double),N,signals);
        //fwrite(channelparity1,sizeof(double),N,signals);
        //fwrite(channelparity2,sizeof(double),N,signals);


		numiter = 
                    decode(channelmesg, channelparity1, channelparity2,N,mesg);

		unsigned numerrors=0;
        unsigned glitches=0;
        unsigned betters=0;


        //bool* comp_msg = new bool[N];
        //fread(comp_msg,sizeof(bool),N,sigs);
        
        
        unsigned loop_number;
        //fread(&loop_number,sizeof(unsigned),1,sigs);
        //if(loop_number>160)
        //    loop_number--;

        for (int x=0;x<N;x++)
        {
            bool temp = channelmesg[x] == 1 ? true : false;

            if (mesg[x] != temp){
                numerrors++;
                /*if (temp != comp_msg[x]){
                    betters++;
                }*/
            }/*else{
                if (temp != comp_msg[x]){
                    glitches++;
                }
            }*/
//            comp_msg[x] = temp;
        }
        //fwrite(comp_msg,sizeof(bool),N,signals);
        //fwrite(L[0],sizeof(double),N,signals);
        //fwrite(L[1],sizeof(double),N,signals);
        //fwrite(&numloops,sizeof(unsigned),1,signals);
        //fwrite(&numerrors,sizeof(unsigned),1,signals);
        //printf("\n%d: errors: %d, glitches: %d, betters: %d",numloops, numerrors, glitches, betters);
      
//        delete comp_msg;

	      	totalerrors += numerrors;
            totalglitches += glitches;
            totalbetters += betters;
	      	totalN += N;
         	totaliter += numiter;
         	numloops++;

	      	char name[30];

		// stopping criterion
		if (totalerrors> 3000)
         		keepgoing = false;
		else if (totalN>maxsim)
         		keepgoing=false;
	}
	while (keepgoing == true); 

	printf ("%5.2lf  %10.7e  %12ld %12ld     %lf     %d %d\n",  dB,
                (totalerrors/(double)totalN), totalerrors, totalN, 
                           (double)totaliter/(double)numloops, N, totaliter);

        fprintf(fp, "%8.4lf %10.7e\n", dB, (totalerrors/(double)totalN));
        fflush(fp);
   }
   fclose(signals);
   //delete interleavearray;
   //delete deinterleavearray;
   //delete mesg;
   //delete parity1;
   //delete parity2;
   //delete channelmesg;
   //delete channelparity1;
   //delete channelparity2;
}




void createencodetable(int states)
{
   bool *boolstate;
   bool *newstate;

if (states == 4) {         // Generators {7,5}
   numstates = 4;
   memory = 2;

   // create arrays used by encode and decode
   output[0] = new bool[numstates];
   output[1] = new bool[numstates];
   fromstate[0] = new unsigned[numstates];
   fromstate[1] = new unsigned[numstates];
   tostate[0] = new unsigned[numstates];
   tostate[1] = new unsigned[numstates];

   boolstate = new bool[memory];
   newstate = new bool[memory];

   for (unsigned input=0;input<2;input++)
      for (unsigned intstate=0;intstate<numstates;intstate++)
      {
         bool boolinput = (input == 0) ? false : true;

         inttobool(intstate,boolstate,memory);

         // calculate output due to the input
         output[input][intstate] = add(boolinput,boolstate[0]);
         output[input][intstate] = add(output[input][intstate],boolstate[1]);

         // calculate new states
         newstate[1] = boolstate[0];
         newstate[0] = add(add(boolinput,boolstate[0]),boolstate[1]);
         // from s' to s
         booltoint (newstate,memory,&tostate[input][intstate]);
         // from s to s'
         fromstate[input][tostate[input][intstate]] = intstate;
      }

}

else if (states == 8) {         // Generators {17,15}
   numstates = 8;
   memory = 3;

   // create arrays used by encode and decode
   output[0] = new bool[numstates];
   output[1] = new bool[numstates];
   fromstate[0] = new unsigned[numstates];
   fromstate[1] = new unsigned[numstates];
   tostate[0] = new unsigned[numstates];
   tostate[1] = new unsigned[numstates];

   boolstate = new bool[memory];
   newstate = new bool[memory];

   for (unsigned input=0;input<2;input++)
      for (unsigned intstate=0;intstate<numstates;intstate++)
      {
         bool boolinput = (input == 0) ? false : true;

         inttobool(intstate,boolstate,memory);

         // calculate output due to the input
         //output[input][intstate] = add(boolinput,boolstate[0]);
         output[input][intstate] = add(boolinput,boolstate[1]);
         output[input][intstate] = add(output[input][intstate],boolstate[2]);
         output[input][intstate] = add(output[input][intstate],boolstate[2]);

         // calculate new states
         newstate[2] = boolstate[1];
         newstate[1] = boolstate[0];
         newstate[0] = add(add(add(boolinput,boolstate[0]),boolstate[1]),
                                                               boolstate[2]);
         // from s' to s
         booltoint (newstate,memory,&tostate[input][intstate]);
         // from s to s'
         fromstate[input][tostate[input][intstate]] = intstate;
      }

}

else if (states == 16) {         // Generators {31,27}
   numstates = 16;
   memory = 4;

   // create arrays used by encode and decode
   output[0] = new bool[numstates];
   output[1] = new bool[numstates];
   fromstate[0] = new unsigned[numstates];
   fromstate[1] = new unsigned[numstates];
   tostate[0] = new unsigned[numstates];
   tostate[1] = new unsigned[numstates];

   boolstate = new bool[memory];
   newstate = new bool[memory];

   for (unsigned input=0;input<2;input++)
      for (unsigned intstate=0;intstate<numstates;intstate++)
      {
         bool boolinput = (input == 0) ? false : true;

         inttobool(intstate,boolstate,memory);

         // calculate output due to the input
         output[input][intstate] = add(boolinput,boolstate[0]);
         output[input][intstate] = add(output[input][intstate],boolstate[3]);
         output[input][intstate] = add(output[input][intstate],boolstate[1]);
         output[input][intstate] = add(output[input][intstate],boolstate[2]);
         output[input][intstate] = add(output[input][intstate],boolstate[3]);

         // calculate new states
         newstate[3] = boolstate[2];
         newstate[2] = boolstate[1];
         newstate[1] = boolstate[0];
         newstate[0] = add(add(boolinput,boolstate[0]),boolstate[3]);
         // from s' to s
         booltoint (newstate,memory,&tostate[input][intstate]);
         // from s to s'
         fromstate[input][tostate[input][intstate]] = intstate;
      }

}

   delete boolstate;
   delete newstate;
}





void inttobool(unsigned state, bool *array, unsigned size)
{
	for (unsigned x=0;x<size;x++)
   {
   	unsigned next = state >> 1;

		if ((next << 1) == state)
      	array[x] = false;
      else
      	array[x] = true;

      state = next;
   }
}




void booltoint(bool *array, unsigned size, unsigned *state)
{
	*state = 0;

   for (int x=0;x<size;x++)
   	if (array[x] == true)
      	(*state) |= (1 << x);
}



unsigned decode (double *mesg, double *parity1, double *parity2, unsigned size, bool *boolmesg)
{
   static double **a[2];
   static double **b[2];
   
   static double **gamma[2][2];
   static double **gammaE[2][2];
   static bool initialized=false;
   static unsigned lastsize=0;
   unsigned returnvalue=ITERATIONS;  

  /* {
	   FILE* ar;
	   ar = fopen("msg.txt","wb");	   
	   fwrite(mesg,sizeof(double),size,ar);			   
	   fwrite(parity1,sizeof(double),size,ar);			   
	   fwrite(parity2,sizeof(double),size,ar);			   
	   fwrite(boolmesg,sizeof(bool),size,ar);
	   fclose(ar);
   }*/

   // minimize new's and delete's to only when needed
   if (size != lastsize && initialized == true)
   {
	   // delete all the arrays and rebuild
	   for (int y=0;y<2;y++)
	   	for (int x=0;x<2;x++)
	      {
	         for (int z=0;z<lastsize;z++)
	         {
	  	      	delete gamma[y][x][z];
	  	      	delete gammaE[y][x][z];
	         }

	      	delete gamma[y][x];
	      	delete gammaE[y][x];
	  		}

	   // create L[encoder #]
	   for (int y=0;y<2;y++)
	      delete L[y];

	   // create alpha[encoder #][k][state]
	   for (int x=0;x<2;x++)
	   {
	   	for (int y=0;y<lastsize;y++)
	      {
	         delete a[x][y];
	         delete b[x][y];
		   }

	      delete a[x];
	      delete b[x];
		}
   }

   if (initialized == false || size != lastsize)
   {
   	initialized = true;
      	lastsize = size;


	// create the arrays dynamically at runtime, delete at end of routine
	// create gamma[encoder #][uk][k][state]
	for (int y=0;y<2;y++)
	 for (int x=0;x<2;x++)
	   {
	    gamma[y][x] = new double*[size];
	      gammaE[y][x] = new double*[size];
	      for (int z=0;z<size;z++)
	      {
	  	 gamma[y][x][z] = new double[numstates];
	         gammaE[y][x][z] = new double[numstates];
	      }
	 }

	// create L[encoder #]
	for (int y=0;y<2;y++)
	   L[y] = new double[size];

	// create alpha[encoder #][k][state]
	for (int x=0;x<2;x++)
	{
	a[x] = new double*[size];
	b[x] = new double*[size];
	// each Yk has 'numstates' values of gamma
	for (int y=0;y<size;y++)
	   {
	   a[x][y] = new double[numstates];
	   b[x][y] = new double[numstates];
	   }
	}
   }

   // initialization of iteration arrays
   for (int x=0;x<numstates;x++)
   	a[0][0][x] = b[0][size-1][x] = a[1][0][x] = (x==0) ? 1.0 : 0.0;

   // initialization of extrinsic information array from decoder 2, 
   // used in decoder 1
   for (int x=0;x<size;x++)
   	L[1][x] = 0.0;

   // 4*Eb/No
   // double Lc = 4.0/No;
   double Lc = 2.0/sigma;

   for (int c=0;c<ITERATIONS;c++)
   {
	// k from 0 to N-1 instead of 1 to N
	for (int k=0;k<size;k++)
	   {
	   // calculate the gamma(s',s);
	   for (int input=0;input<2;input++)
		{
		double uk = (input == 0) ? -1.0 : 1.0;

	      	for (int s=0;s<numstates;s++)
	         {
	         double xk = (output[input][s] == 0) ? -1.0 : 1.0;

		 gammaE[0][input][k][s]=exp(0.5*Lc*parity1[k]*xk);
		 gamma[0][input][k][s]
                      = exp(0.5*uk*(L[1][k]+Lc*mesg[k]))*gammaE[0][input][k][s];
	         }
	      }
	   }

	//debugging loop
	/*{
		FILE* ar;
		ar = fopen("gama0.txt","wb");
		for (int k=0;k<size;k++){			
			for (int input=0;input<2;input++){
				for (int s=0;s<numstates;s++){
					fwrite((void*)&gamma[0][input][k][s],sizeof(double),1,ar);
				}
			}
		}
		fclose(ar);
	}*/

	   // calculate the alpha terms
	   // from 1 to N-1, 0 is precalculated, N is never used
	   for (int k=1;k<size;k++)
	   {
	      double temp=0;
	      // calculate denominator
	      for (int state=0;state<numstates;state++)
	         temp += a[0][k-1][fromstate[0][state]]
                                       * gamma[0][0][k-1][fromstate[0][state]] 
                       + a[0][k-1][fromstate[1][state]]
                                       * gamma[0][1][k-1][fromstate[1][state]];
	      for (int state=0;state<numstates;state++)
		a[0][k][state] = ( a[0][k-1][fromstate[0][state]]
                                       * gamma[0][0][k-1][fromstate[0][state]] 
                                 + a[0][k-1][fromstate[1][state]]
                                       * gamma[0][1][k-1][fromstate[1][state]] )
                                 / temp;
	   }

	   /*{
		   FILE* ar;
		   ar = fopen("alpha0.txt","wb");
		   for (int k=0;k<size;k++){		
			   for (int s=0;s<numstates;s++){
				   fwrite((void*)&a[0][k][s],sizeof(double),1,ar);
			   }
		   }
		   fclose(ar);
	   }*/

	   // calculate the beta terms
	   // from N-1 to 1
	   for (int k=size-1;k>=1;k--)
	   {
	      double temp=0;
	      // calculate denominator
	      for (int state=0;state<numstates;state++)
		temp += a[0][k][fromstate[0][state]]
                                       * gamma[0][0][k][fromstate[0][state]] 
                      + a[0][k][fromstate[1][state]]
                                       * gamma[0][1][k][fromstate[1][state]];
	      for (int state=0;state<numstates;state++)
	      	b[0][k-1][state] = ( b[0][k][tostate[0][state]]
                                       * gamma[0][0][k][state] 
                                   + b[0][k][tostate[1][state]]
                                       * gamma[0][1][k][state]  ) / temp;
	   }

	   /*{
		   FILE* ar;
		   ar = fopen("beta0.txt","wb");
		   for (int k=0;k<size;k++){		
			   for (int s=0;s<numstates;s++){
				   fwrite((void*)&b[0][k][s],sizeof(double),1,ar);
			   }
		   }
		   fclose(ar);
	   }*/


	   // compute the extrinsic LLRs
	   for (int k=0;k<size;k++)
	      {
	      double min=0;

	      // find the minimum product of alpha, gamma, beta
	      for (int u=0;u<2;u++)
	         for (int state=0;state<numstates;state++)
	         {
	            double temp = a[0][k][state]
                       * gammaE[0][u][k][state] * b[0][k][tostate[u][state]];
	            if ((temp < min && temp != 0)|| min == 0)
	         	min = temp;
	         }

	      // if all else fails, make min real small
	      if (min == 0 || min > 1)
	         min = 1E-100;

	      double topbottom[2];
	      for (int u=0;u<2;u++)
	         {
		   topbottom[u]=0.0;
	   	   for(int state=0;state<numstates;state++)
	        	topbottom[u] += 
                            (a[0][k][state] * gammaE[0][u][k][state]
                                            * b[0][k][tostate[u][state]] );
	         }

	      if (topbottom[0] == 0)
	         topbottom[0] = min;
	      else if (topbottom[1] == 0)
	         topbottom[1] = min;

	      L[0][k] = (log(topbottom[1]/topbottom[0]));
         }

	   /*{
		   FILE* ar;
		   ar = fopen("llr0.txt","wb");
		   for (int k=0;k<size;k++){					   
			   fwrite((void*)&L[0][k],sizeof(double),1,ar);			   
		   }
		   fclose(ar);
	   }*/

	interleavedouble(L[0],size);
	// remember to deinterleave for next iteration
	interleavedouble(mesg,size);

	// start decoder 2
	// code almost same as decoder 1, 
        // could combine code into one but too lazy
	for (int k=0;k<size;k++)
	   {
	   // calculate the gamma(s',s);
	   for (int input=0;input<2;input++)
		{
		double uk = (input == 0) ? -1.0 : 1.0;

	      	for (int s=0;s<numstates;s++)
	            {
	            double xk = (output[input][s] == 0) ? -1.0 : 1.0;

		    gammaE[1][input][k][s] = exp(0.5*Lc*parity2[k]*xk);
		    gamma[1][input][k][s] = exp(0.5*uk*(L[0][k]+Lc*mesg[k]))
                                                     * gammaE[1][input][k][s];
	            }
	         }
	   }

		//debugging loop
	/*{
		FILE* ar;
		ar = fopen("gama1.txt","wb");
		for (int k=0;k<size;k++){			
			for (int input=0;input<2;input++){
				for (int s=0;s<numstates;s++){
					fwrite((void*)&gamma[1][input][k][s],sizeof(double),1,ar);
				}
			}
		}
		fclose(ar);
	}*/

	// calculate the alpha terms
	for (int k=1;k<size;k++)
	   {
		double temp=0;

	   	// calculate denominator
	      	for (int state=0;state<numstates;state++)
		     	temp += a[1][k-1][fromstate[0][state]]
                                     * gamma[1][0][k-1][fromstate[0][state]] 
                              + a[1][k-1][fromstate[1][state]]
                                     * gamma[1][1][k-1][fromstate[1][state]];

	     	for (int state=0;state<numstates;state++)
			a[1][k][state] = (a[1][k-1][fromstate[0][state]] 
                                       * gamma[1][0][k-1][fromstate[0][state]]
                                        + a[1][k-1][fromstate[1][state]]
                                       * gamma[1][1][k-1][fromstate[1][state]] )
                                                   / temp;
	   }

	 /*{
		   FILE* ar;
		   ar = fopen("alpha1.txt","wb");
		   for (int k=0;k<size;k++){		
			   for (int s=0;s<numstates;s++){
				   fwrite((void*)&a[1][k][s],sizeof(double),1,ar);
			   }
		   }
		   fclose(ar);
	   }*/


	// in the first iteration, set b[1][N-1] = a[1][N-1] for decoder 2
	// this decoder can't be terminated to state 0 because of the
	// interleaver. The performance loss is supposedly neglible.
      	if (c==0)
	   {
		double temp=0;

		// calculate denominator
		for (int state=0;state<numstates;state++)
			temp += 
                          a[1][size-1][fromstate[0][state]]
                                    * gamma[1][0][size-1][fromstate[0][state]] 
                          + a[1][size-1][fromstate[1][state]]
                                    * gamma[1][1][size-1][fromstate[1][state]];

		for (int state=0;state<numstates;state++)
			b[1][size-1][state] = 
                          (a[1][size-1][fromstate[0][state]]
                                    * gamma[1][0][size-1][fromstate[0][state]] 
                         + a[1][size-1][fromstate[1][state]]
                                    * gamma[1][1][size-1][fromstate[1][state]] )
                                                    / temp;
	   }

	for (int k=size-1;k>=1;k--)
	   {
	   	double temp=0;

	   	// calculate denominator
	      	for (int state=0;state<numstates;state++)
		     	temp += a[1][k][fromstate[0][state]]
                                 * gamma[1][0][k][fromstate[0][state]] 
                              + a[1][k][fromstate[1][state]]
                                 * gamma[1][1][k][fromstate[1][state]];

	   	for (int state=0;state<numstates;state++)
	      		b[1][k-1][state] = 
                         (b[1][k][tostate[0][state]]*gamma[1][0][k][state] 
                           + b[1][k][tostate[1][state]]*gamma[1][1][k][state])
                                                   / temp;
	   }
	   /*{
		   FILE* ar;
		   ar = fopen("beta1.txt","wb");
		   for (int k=0;k<size;k++){		
			   for (int s=0;s<numstates;s++){
				   fwrite((void*)&b[1][k][s],sizeof(double),1,ar);
			   }
		   }
		   fclose(ar);
	   }*/

	for (int k=0;k<size;k++)
	   {
	   	double min = 0;

	      	// find the minimum product of alpha, gamma, beta
	      	for (int u=0;u<2;u++)
	      	for (int state=0;state<numstates;state++)
	         {
	         	double temp=a[1][k][state]
                                         *gammaE[1][u][k][state]
                                                   *b[1][k][tostate[u][state]];

	            	if ((temp < min && temp != 0)|| min == 0)
	            	   min = temp;
	         }
	      	// if all else fails, make min real small
	      	if (min == 0 || min > 1)
	      	   min = 1E-100;

		double topbottom[2];

		for (int u=0;u<2;u++)
		{
		   	topbottom[u]=0.0;

	   	   	for(int state=0;state<numstates;state++)
			   topbottom[u] += (a[1][k][state]
                                            * gammaE[1][u][k][state]
                                            * b[1][k][tostate[u][state]]);
		}

	      	if (topbottom[0] == 0)
	      	   topbottom[0] = min;
	        else if (topbottom[1] == 0)
	      	   topbottom[1] = min;

	      	L[1][k] = (log(topbottom[1]/topbottom[0]));
	   }

		/*{
		   FILE* ar;
		   ar = fopen("llr1.txt","wb");
		   for (int k=0;k<size;k++){					   
			   fwrite((void*)&L[1][k],sizeof(double),1,ar);			   
		   }
		   fclose(ar);
	   }*/

	deinterleavedouble(mesg,size);
	deinterleavedouble(L[1],size);
	// get L[0] back to normal after decoder 2
	deinterleavedouble(L[0],size);

      	bool temp=true;
		int e=0;
		/*bool* temp_resp = new bool[size];
		double* soft_bit = new double[size];*/

		for (int k=0;k<size;k++)
		{
			if (boolmesg[k] != (( Lc*mesg[k]+L[0][k]+L[1][k]>0.0 ) ? true:false ))
			{
				temp = false;	
				e++;
			}
//			temp_resp[k] = boolmesg[k] != (( Lc*mesg[k]+L[0][k]+L[1][k]>0.0 ) ? true:false );
//			soft_bit[k] = ( Lc*mesg[k]+L[0][k]+L[1][k] );
		}

      	// we can quit prematurely since it has been decoded
      	if (temp==true)
      	{
      		returnvalue = c+1;
		c = ITERATIONS;
      	}
   }
   // end decoder 2

   // make decisions
   for (int k=0;k<size;k++)
	if ((Lc*mesg[k] + L[0][k] + L[1][k]) > 0)
      		mesg[k] = 1.0;
      	else
      		mesg[k] = -1.0;

   return returnvalue;
}




void deinterleavedouble(double *mesg, unsigned size)
{
	double *temp;

   temp = new double[size];

   for (int x=0;x<size;x++)
   	temp[x] = mesg[x];

	for (int x=0;x<size;x++)
		mesg[deinterleavearray[x]] = temp[x];

   delete temp;
}




void interleavedouble(double *mesg, unsigned size)
{
	double *temp;

   temp = new double[size];

   for (int x=0;x<size;x++)
   	temp[x] = mesg[x];

   for (int x=0;x<size;x++)
   	mesg[interleavearray[x]] = temp[x];
   delete temp;
}




void interleave(bool *mesg,unsigned size)
{
	bool *temp;

   temp = new bool[size];

   for (int x=0;x<size;x++)
   	temp[x] = mesg[x];

   for (int x=0;x<size;x++)
   	mesg[interleavearray[x]] = temp[x];

   delete temp;
}




void deinterleave(bool *mesg,unsigned size)
{
	bool *temp;

   temp = new bool[size];

   for (int x=0;x<size;x++)
   	temp[x] = mesg[x];

	for (int x=0;x<size;x++)
		mesg[deinterleavearray[x]] = temp[x];

   delete temp;
}








void encode(bool *mesg,bool *parity,unsigned size, bool force)
{

	unsigned state=0;

  	for (int x=0;x<size;x++)
	{
		// force the encoder to zero state at the end
     	if (x>=size-memory && force)
      {
			if (tostate[0][state]&1)
         	mesg[x] = true;
         else
         	mesg[x] = false;
      }

		// can't assume the bool type has an intrinsic value of 0 or 1
		// may differ from platform to platform
      int uk = mesg[x] ? 1 : 0;

		// calculate output due to new mesg bit
      parity[x] = output[uk][state];
      // calculate the new state
		state = tostate[uk][state];
	}
}




bool add(bool a, bool b)
{
   return a==b ? false : true;
}

//double gaussian(double sigma)
//{
//  double rndm, u1, u2, s, noise;
//    do {
//        rndm = (double)(random())/MAXRAND;
//        u1 = rndm * 2 - 1.0;
//        rndm = (double)(random())/MAXRAND;
//        u2 = rndm * 2 - 1.0;
//        s = u1 * u1 + u2 * u2;
//      } while( s >= 1 );
//    noise = u1 * sqrt( (-2.0*log(s))/s );
//  return(noise*sigma);
//}

double 
gaussian(double sigma)
{
  double rndm, u1, u2, s, noise;
    do {
        rndm = (double)(rand())/RAND_MAX;
        u1 = rndm * 2 - 1.0;
        rndm = (double)(rand())/RAND_MAX;
        u2 = rndm * 2 - 1.0;
        s = u1 * u1 + u2 * u2;
      } while( s >= 1 );
    noise = u1 * sqrt( (-2.0*log(s))/s );
  return(noise*sigma);
}




